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1 .  INTRODUCTION 

Computational  fluid  dynamics  (CFD)  is  a  tool  which 
predicts  the  gas  dynamics  of  blast  problems  of  interest  to  the 
Army  by  solving  a  set  of  mathematical  equations  with  a 
high-speed  digital  computer.  The  governing  equations  for  the 
blast  problem  presented  here  are  the  two-dimensional  unsteady 
Euler  equations.  The  computations  were  performed  on  a  Cray 
XMP/48  supercomputer  by  discretizing  the  Euler  equations  with 
an  upwind,  Total  Variation  Diminishing  (TVD) ,  finite-volume, 
implicit  scheme.  In  a  paper  by  Molvik.(l)*  the  scheme  was 
presented  in  detail  and  proved  to  be  well  suited  for  blast 
wave  calculations.  The  scheme  is  discussed  in  the 
computational  algorithm  section.  The  algorithm  is  used  here 
to  provide  gas  dynamic  information  for  a  candidate  Large-Scale 
Blast  and  Thermal  Simulator  (LB/TS)  concept. 

The  Army  has  a  growing  need  for  nuclear  blast  and  thermal 
survivability  testing  of  tactical  equipment.  In  order  to  meet 
this  need,  the  Army  is  conducting  research  into  the  design  and 
operation  of  a  Large-Scale  Blast  and  Thermal  Simulator  (LB/TS), 
essentially  a  large  multi-driver  shock  tube  with  thermal 
capabilities.  The  LB/TS  deslgn(2)  currently  consists  of  a 
number  of  driver  tubes  releasing  compressed  gas  through  a 
series  of  converging-diverging  nozzles  into  a  large  expansion 
tunnel.  Figure  1.  The  compressed  gas  forms  a  shock  which 
travels  down  the  expansion  tunnel  and  produces  the  blast 
simulation.  The  thermal  simulation  is  accomplished  by  igniting 
a  Thermal  Radiation  Source  (TRS)  based  on  aluminum /oxygen 
combustion  Just  before  the  arrival  of  the  blast  wave  at  the 
test  section.  The  expansion  tunnel  is  physically  large  enough 
to  accommodate  the  testing  of  full-sized  tactical  vehicles  such 
as  tanks  and  helicopters  at  low  blockage  of  the  test  section. 

The  LB/TS  can  be  modeled  in  a  2-D  axisymmetric  sense  by 
combining  the  drivers  and  nozzles  to  form  one  equivalent  driver 
as  shown  in  Figure  3.  This  simplified  model  of  the  LB/TS  was 
actually  built  for  experimental  testing  at  Aberdeen  Proving 


*  The  numbers  inside  (  )  denote  the  reference  number. 
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Ground.  The  thermal  radiation  source  was  not  experimentally  or 
computationally  modeled.  This  shock  tube  is  l/5?th  of  the  scale 
of  the  proposed  LB/TS.  The  experimental  data  obtained  in  this 
tube  was  used  for  LB/TS  design  studies  and  BLAST2D  code 
validation.  All  of  the  experimental  and  computational  data 
presented  in  this  paper  was  from  this  shock  tube 
configuration.  Blast  waveforms  were  produced  with  peak 
overpressures  ranging  from  approximately  5  psi  to  35  psi. 
Heating  of  the  driver  gas  was  performed  for  some  of  the  high 
pressure  cases  to  reduce  the  driver  pressure  required  to  obtain 
a  given  shock  overpressure,  and  to  alleviate  the  temperature 
(density)  discontinuity  at  the  contact  surface  between  the 
expanded  driver  gas  and  shocked  expansion  section  gas. 

Currently,  one-dimensional  calculations  have  been 
performed  for  the  1/57  scale  LBS  with  useful  results . (3 , 4) 
However,  the  one-dimensional  calculations  have  had  limited 
success  for  accurately  predicting  the  flow  through  the 
diverging  portion  of  the  LBS  design  because  the  flow  in  this 
region  is  multi-dimensional.  The  flow  is  multi-dimensional 
due  to  the  rapid  and  large  area  change  that  exists  in  the 
diverging  nozzle.  The  remainder  of  this  paper  presents  the 
upwind,  TVD,  finite-volume,  Implicit  scheme  in  the  BRL  BLAST2D 
code  and  results  vhich  capture  and  reveal  the  nature  of  the 
flow  physics  in  the  1/57  scale  LBS. 


2 .  GOVERNING  EQUATIONS 

The  governing  equation  is  the  Euler  equation  written  in 
integral  form: 

d  C  f  ,  _ 


f  Qdl/  +  f  n  •  FdS  = 

•'ll  j  c 


The  integral  form  of  the  Euler  equation  can  be  rewritten 
for  a  two-dimensional  generalized  cell  volume  (Figure  3)  as  : 
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where ; 


Q  = 


E 


'  pU  \ 
pUu  +  Y^p  \ 

pUv  -  x^p  1 

(e+p)U  J 


F  = 


pV 

pVu  -  y?p 
pVv  +  x^p 
(e+p)V 


(3) 


This  set  of  four  integral  equations  represents  the 
conservation  of  mass,  momentum  in  x  and  y  directions,  and 
energy,  per  unit  volume  where  p  is  the  density,  p  is  the 
pressure,  u  and  v  are  the  velocities  in  the  x 
(longitudinal)  and  y  (height)  directions  respectively, 
and  e  is  the  total  internal  energy  per  unit  volume: 

e  =  )  +  V2  p(u2  +  v2) 


The  volume  fluxes  are  defined  as: 


U  =  y  u  -  x  v  (5 

n  n 

V  =  yru  +  xrv 

5  *  (6 


For  the  two  dimensional  cell  shown  in  Figure  3,  the 
integration  of  flux  over  the  surface  in  Equation  1  has  been 
replaced  in  Equation  2  by  an  integral  over  each  face  of  the 
cell.  The  P-direction  is  taken  as  the  body  normal  and  the 
C-direction  is  tangential  to  the  surface  of  the  body.  The  cell 
volume  and  walls  are  assumed  to  be  fixed  in  time.  The  metrics 


are  the  vector  elements  of  the  cell  walls  and  V 


is  the  volume  of  the  grid  cell. 


The  physical,  independent  variables  (x,y,t)  were 
transformed  into  a  uniformly  spaced  computational  grid  (s.n.-x) 
by  a  general  transformation  of  the  form: 


.*  •>»  f+m •  ‘ 


T  =  t 

C  =  ;(t.x,y) 
n  =  n(t,x,y) 


The  transformations  were  chosen  so  that  the  grid  spacing  in 
the  computational  space  is  uniform  and  of  unit  length,  ac  -  l , 
An=  1.  Thus,  the  uniform  equi-spaced  mesh  in  psi  and  eta 
allows  the  use  of  unweighted  differencing  schemes.  As  a 
result ,  the  computational  code  can  be  applied  to  a  variety  of 
physical  geometries  and  grids. 

If  an  average  flux  is  defined  on  the  cell  faces,  and  ac 
and  An  are  taken  as  unity,  the  integral  form  of  the  Euler 
equation,  Equation  2  can  be  rewritten  in  finite  volume  form 


Q0!  -  q"  . 
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where  the  indices  i  and  J  correspond  to  the  C  and  n  directions 
respectively  in  the  computational  mesh  as  shown  in  Figure  3 . 
The  vectors  E  and  F  are  the  convective  numerical  fluxes  in 
computational  space  (c.n.t)  consistent  with  the  physical 
fluxes  E  and  F  in  (x,y,t).  The  vector  Q  consists  of  the  cell 
averaged  dependent  variables.  The  integration  scheme  is  fully 
implicit  if  m=n+l  and  is  explicit  if  m=n.  The  variables  have 
been  nondimenslonalized  as  follows ; 
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where  L=l,  c=sound  speed,  and  subscripts  1  and  4  represent 
conditions  in  the  driven  and  driver  sections  respectively. 


3.  THE  COMPUTATIONAL  ALGORITHM 


3.1  Introduction.  Discretization  of  the  governing 
equations  into  an  upwind,  TVD,  finite-volume,  implicit  scheme 
produces  an  algorithm  that  is  well  suited  for  blast  wave 
calculations. Cl)  Upwind  flux  difference  splitting  with  TVD 
achieves  second-order  accuracy  without  introducing  spurious 
oscillations  near  discontinuities.  Strong  gradients  and 
complex  flow  fields  are  resolved  accurately.  Older  techniques 
used  central  differencing  schemes  with  arbitrary  smoothing 
parameters  which  could  not  be  relied  on  to  capture  strong 
gradients  (i.e. ,  pressure  ratio  across  the  shock  >  10.0) 
accurately . (5)  The  advantages  of  the  central  differencing 
techniques  were  programming  simplicity  and  adequate  resolution 
for  weak  gradient  problems.  However,  for  the  complex  flow 
fields  and  strong  gradients  typical  of  blast  problems  upwind 
differencing  with  TVD  provides  better  resolution.  The 
disadvantages  of  upwind  differencing  with  TVD  are  long 
computing  times  caused  by  an  increase  in  the  number  of 
arithmetic  operations  per  integration  step  and  loss  of 
programming  simplicity.  The  BLAST2D  code  results  shown  in  this 
paper  were  generated  on  a  Cray  XMP/48  and  typically  took  six 
to  seven  hours  of  cpu  time. 


Conservative  schemes  capture  shocks  and  other 
discontinuities  automatically.  The  finite  volume  philosophy 
ensures  conservation  at  interior  and  boundary  points.  The 
scheme  is  made  implicit  by  linearizing  only  the  first-order 


contribution  and  by  employing  a  Newton  iteration  of  the  type 
described  by  Rai(5)  to  eliminate  any  approximations  made.  The 
Implicit  version  of  the  scheme  requires  more  computations  per 
integration  step  than  the  explicit  version,  but  permits  larger 
time  steps  which  overall  reduces  computational  expense. 

The  next  section  presents  the  first-order  accurate  upwind 
schemed)  which  is  the  basic  building  block  of  the 
computational  algorithm.  Subsequently,  the  first-order  scheme 
is  expanded  to  second  order  accuracy  with  the  addition  of 
second-order  terms  and  TVD  concepts.  Development  of  the 
implicit  version  of  the  algorithm  and  the  Newton  iterative 
procedure  used  is  presented.  Finally,  boundary  conditions  are 
discussed  briefly. 

3.2  First-Order  Scheme.  The  first-order  scheme  is 
based  upon  Roe's  approximate  Riemann  solver  (6,?)  coupled  with 
upwind  flux  difference  splitting.  First,  approximate  Riemann 
solvers  are  discussed.  Then,  the  information  supplied  by  the 
Riemann  solver  is  used  with  upwind  flux  difference  splitting 
concepts  to  provide  the  first-order  convective  fluxes  E  and  F 
in  the  finite  volume  form  of  the  Euler  equation,  Equation  8. 

Riemann  problems  are  incorporated  into  the  numerical 
solution  by  considering  the  dependent  variables  at  cell 
centers  for  each  cell  in  turn,  as  pairs  of  states  defining  a 
sequence  of  Riemann  problems,  Figure  4.  The  Riemann  problem 
for  the  {.  direction  in  Figure  4  is  :  given  two  states 
(pl.ul.pl)  and  (p4,u4,p4)  determine  the  combination  of  shocks, 
contact  discontinuities,  and  expansions  which  result  in  these 
end  states,  that  is,  determine  (p2,u2,p2)  and  (p3,u3,p3).  To 
obtain  a  solution,  exact  Riemann  solvers  require  an  iterative 
procedure  which  is  computationally  expensive  when  performed 
for  a  large  number  of  cells  and  time  steps.  The  expense  of 
producing  an  exact  solution  to  the  Riemann  problem  is 
justified  only  if  the  information  made  available  could  be  put 
to  some  sophisticated  use.  The  approximate  Riemann  solvers  are 
considerably  less  expensive  because  the  Riemann  problem  is 
solved  with  a  direct  non-iterative  method  which  is  about  as 
time  consuming  as  one  cycle  of  the  iterative  procedures. 
Comparisons  of  the  solutions  from  the  exact  vs.  approximate 
Riemann  solvers  reveal  slight  differences.  Other  approximate 
Riemann  solvers  could  have  been  used,  but  Roe's  method  is 
the  approach  recommended  by  Chakravarthy  when  computational 


res* 


'«».  ■  .V  »  »  -  -  •  o  ■>  < 


S’  %' 


efficiency  is  important . (6) 

References  6  and  7  outline  in  detail  Roe's  method 
for  determining  the  intermediate  states  of  the  Riemann 
problem.  In  general,  the  solution  consists  of  four 
constant  property  states  separated  by  three  waves,  Figure 
4.  Once  the  dependent  variables  are  obtained  at  the 
intermediate  states,  the  flux  at  the  cell  interface  is 
calculated  by  determining  the  flux  change  across  the 
waves.  The  flux  change  associated  with  the  waves 
traveling  in  the  positive  £  direction  is  given  the 
symbol  AE+and  that  in  the  negative  direction  is 
represented  by  ae~.  The  waves  carry  information  from 
the  "upwind"  direction  to  the  cell  center,  thus  the  notion 
of  upwind  differencing.  The  flux  remaining  at  the  inter¬ 
face  for  all  time  associated  with  this  Riemann  problem  must 
then  be  represented  by  either  of  the  following  equations: 


Ei+>S  =  Ei  +  4EH>5 


Ei+!i  =  Ei  +  t  ■  4Ei+is 

or,  by  averaging  the  two  equations  above, 

Ei+!s  =  1/2  (Ei  +  ehi  +  4EUS  -  <*s) 


(10) 

(11) 


(12) 


Let  A  ,  R,  and  R  denote  the  eigenvalue  matrix,  and  the 
right  and  left  eigenvector  matrices  respectively,  evaluated  at 
the  cell  interface . (7)  The  flux  difference  across  the  positive 
and  negative  velocity  waves  can  be  calculated.  They  are: 


(13) 


However,  the  dependent  variables  are  not  defined  at  the  cell 
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interfaces  where  these  matrices  must  be  evaluated.  Roe (7)  has 
developed  a  special  averaging  process  to  calculate  the 
dependent  variables  on  the  cell  interface  and  satisfy  the 
following  relation. 


*,♦,  -E,  ■  MS!  (*i*i  -  *1) 


(15) 


The  superscript  Roe  denotes  a  Roe  averaged  quantity.  By 
satisfying  the  relations  above,  the  shock  capturing 
capabilities  of  the  algorithm  are  retained  and  correct  wave 
speeds  are  assured.  Roe's  averaging  of  the  dependent 
variables  proceeds  as  follows : 

uiV^7+  _  viV^7  +  viV^T 

i+1*  yfti 


hi+1!  = 


_  hiVpi +  hi+iVpuT 


y/°u- 


■UH  -  j(hHb  -  1/2  (vh  +  »U)) 


where  the  total  enthalpy  per  unit  mass  is 


h  =  (e  +  p)/p 


(16) 


(17) 


The  first-order  flux  on  the  j+1/2  interface  can  be  obtained  in 
a  similar  manner  by  replacing  x.  with  -  xp  and  y ^  with  -yr  - 

3.3  Second-Order  Scheme.  A  second-order  convective 
flux  can  be  produced  by  adding  a  correction  term  to  the 
first-order  flux.  However,  in  order  to  avoid  spurious 
oscillations,  the  correction  term  must  fulfill  the  criteria 
for  the  algorithm  to  be  TVD.  TVD  schemes  achieve  second-order 
accuracy  without  introducing  spurious  oscillations  near 
discontinuities  by  employing  a  feedback  mechanism  -"smart 
numerical  dissipation"-  wherein  fluxes  are  compared  at 
neighboring  control  volumes.  In  regions  of  little  change  no 
numerical  dissipation  is  added  to  the  second  order  correction 
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terms,  while  in  regions  of  large  change,  numerical  dissipation 
is  added  to  ensure  stability. 

During  this  process  no  new  extrema  are  created  by  the 
numerical  dissipation.  TVD  data  preserve  monotonicity;  a)  no 
new  extrema  must  be  created  and  b)  the  absolute  value  of  any 
extrema  must  not  increase.  TVD  schemes  yield  oscillation-free 
solutions  by  modifying  flux  differences  to  meet  the  above 
criteria.  Reference  7  outlines  a  class  of  explicit  flux 
limiting  schemes  that  fulfill  this  criteria.  The  second-order 
flux  for  the  fully  upwind  scheme  can  be  written  as(l): 


♦  1/2 


[afU  -  4iW] 


(18) 


If  the  following  definitions  are  made, 
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(19) 
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(20) 

then  the  limited  values  of  the  flux  can  be  written  as 


^i+%  Ri+Js  Aai+^  ’  =  Ri+4  Aai+!j 


(21) 


Note  that  the  characteristic  fluxes  are  limited,  not  the 
fluxes  given  in  Equation  3.  The  symbols  ~  and  «  shown  over 
the  to  denote  flux-limited  values  of  Ao  and  are  computed  as 
follows : 


AoT+v,  =  minmod  ["  AoT+l,  ,  BAo^J 
Aa|+is  =  minmod  ["  Ao*+Jj  ,  BAat+3/2J 


(22) 

(23) 


Where  the  "minmod"  is  defined  as 

minmod  [x,y]  =  sign  (*)  *  max  (0,  minjj*|,  ysign  (x)}]  (24) 

and  beta  is  a  compression  parameter  that  is  restricted  to  fall 
in  the  range 
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1  <  6  <  2 


(25) 


Equation  8  can  be  rewritten  with  the  first-order 
convective  fluxes,  E  and  P,  replaced  with  the  second-order 
f luxes : 


qH+1  qA  r2nd  r2nd 

j nLl  i.J  +  1 I  hhA 


hi+k 


p2nd 

i 


=  0  (26) 


The  above  discussion  describes  the  explicit  second-order 
accurate  in  space  scheme.  Second-order  accuracy  in  time  is 
achieved  by  simply  replacing  the  first-order,  backward 
derivative  of  the  time-dependent  variables  with  a  second-order 
backward  difference. 

3.4  Implicit  Scheme.  For  a  fully  implicit  scheme, 
the  fluxes  must  be  evaluated  at  the  n+1  time  level.  In  order 
to  calculate  a  flux  at  the  n+1  time  level  the  flux  must  be 
linearized  with  respect  to  time,  t.  The  first-order  numerical 
flux  on  the  i+1/2  cell  interface  evaluated  at  the  n+1  time 
level  is  represented  as: 


Fn+1  .  I  cn+l 
i+H  '2  -i+1 


r1  +  <a"  -  a+c2<!  -  of) 


(27) 


An  approximate  linearization  of  this  interface  flux  may  be 
achieved  by  freezing  the  coefficient  matrix  (A”  -  A+)  at  time 
level  n  and  linearizing  the  remaining  terms.  Numerical 
experiments  have  shown  that  such  an  approximation  is 
acceptable.  The  linearized  numerical  flux  is  then  written  as: 


rn+1  =  1  An 
i+h;  2  i 


A?+l  +  (A'  -  A+>?*  +  2  A?  ■  -  A+>U  -><>1  *  E? 


+H  (28) 


=  (AK)?  AQj+)  +  (flL)"  iQ(  *  E 1 


where , 


AQj  *  Q^+1  - 


Using  a  similar  type  of  linearization  for  the  body  normal  flux, 
F,  as  for  the  streamwise  flux,  E,  the  first-order,  implicit 
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numerioal  algorithm  is  written  as: 
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At  J/rn  ;n  ,  ,  /in  ;n  s 

V-  i  L  Ei+i5J  ~  Fi»j+%  "  Fi,j-% 

*  5  J 

To  avoid  the  expense  of  inverting  a  large  sparse  matrix,  an 
approximate  faotorization  is  done  to  break  the  banded  matrix 
into  two  tridiagonal  matrioes .  This  is  written  in  two  steps 
with  the  asterisked  *  variables  denoting  an  intermediate  step 
as : 

* 

-  -  ‘aLi?v5*-i.j]  ‘ RHS(32)  (30> 


>(d  +  v:[(BR,j^id*i +  ubL,U  -  ,BR>3V45>d 
^"^id-ij-^d 


(31) 


3-5  Iterative  Soheme.  In  order  to  eliminate  the 
linearization  and  approximate  factorization  errors  that  might 
occur,  a  Newton  iteration  technique  is  employed.  The 
iteration  takes  the  form: 

i5*.j +  v\-  [(flR)W5iVu  +  1(flL)U  -  (AR)?-H)SB*,j  - 

i  >  j  _  J  ’ 


Kd  -  -  vf:  (^i-^d)t(f5d«,-f?d--J 

1  9  J  L 


(32) 
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where  aQ  is  now  defined  as  the  iterative  change  in  the 
cell  averaged  dependent  variables, (  Q?+l  -  Q?  •  )rather 
than  the  time  change  and  p  denotes  the’ iteration  number. 
Ideally,  all  linearization  errors  are  completely 
eliminated  when  the  residual  of  equation  31  is  driven  to 
zero.  However,  in  practice,  convergence  is  defined  short 
of  this  with  minimal  loss  in  accuracy,  in  order  to  reduce 
the  number  of  iterations  and  hence  the  expense  of  the 
calculation. 


3.6  Axlsymmetrlc  Source  Term.  The  governing  equations 
and  numerical  algorithm  described  so  far  are  appropriate  for 
2-D  planar  problems.  It  is  a  simple  matter  to  Include 
axisymmetric  effects  by  adding  the  source  terms  denoted  by  H 
and  I  to  equation  32  where: 
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Thus,  equations  32  and  33  become: 
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3.7  Boundary  Conditions.  The  invlsoid  boundary 
conditions  are  obtained  by  specifying  an  appropriate  flux  on 
the  walls  of  the  shock  tube  and  at  the  symmetry  boundary.  Only 
half  of  the  symmetrical  shock  tube  is  actually  modeled  to  save 
computational  expense.  Then,  the  results  are  graphically 
reflected.  At  the  walls  and  symmetry  boundary,  the  normal 
component  of  velocity  is  zero,  the  tangential  component  of 
velocity  is  nonzero.  The  flux  on  these  surfaces  can  then 
be  represented  as: 


F  =  -V 


(38) 


Only  a  value  of  pressure  need  be  evaluated  at  the 
surface.  As  a  first  approximation,  one  might  consider  using 
the  pressure  of  the  cell  directly  above  the  surface.  This 
translates  into  a  zero-order  approximation.  However  a 
first-order  approximation  of  the  surface  flux  can  be  made  if  a 
Riemann  problem  is  set  up  on  the  surface.  This  is  consistent 
with  the  interior  scheme  and  would  seem  like  the  reasonable 
approach.  The  first-order  Riemann  solver  is  used  between 
the  first  cell  off  the  surface  and  a  reflected  cell.  If 
the  subscripts  1  and  -1  denote  the  first  cell  off  the 
surface  and  the  reflected  cell  respectively,  the  surface  flux 
can  then  be  written  as: 


K  =  1/2 


_F1  +  F-!  +  (fl‘  '  4  («,  -  ?,)" 


(39) 
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The  dependent  variables  of  the  reflected  cell  are  calculated 
using  the  following  relations. 


P-1  *  P1 


-i  =  h 


“-1  =  (*S2  -  yS2)Ul  +  2V?Vl]/[xJ2  +  yS2l  (40) 

v-i =  (y2  -  Xs2)', +  2v5uij/  V  +  yt2] 

The  metrics  in  the  above  equations  are  those  of  the  cell 
interface  on  the  surface.  A  second  order  flux  can  be  obtained 
by  reflecting  even  another  set  of  dependent  variables  with  a 
subscript  of  -2. 

4.  GEOMETRY,  GRID,  AND  INITIAL  CONDITIONS  FOR  LBS  COMPUTATION 


The  geometry  of  the  single-driver,  1/57  scale  model  LBS 
design  concept  is  shown  in  Figure  2.  The  corners  and  sharp 
edges  were  smoothed  out  in  the  computational  shock  tube  to 
simplify  grid  generation  and  to  facilitate  the  use  of  an 
elliptic  grid  generator  for  the  driver  and  converging- 
diverging  nozzle  portion  of  the  tube.  The  elliptic 
grid  generator  produced  grid  cells  that  varied  smoothly  in 
regions  of  rapid  change  such  as  the  converging-diverging 
nozzle.  The  grid  generated  for  this  configuration  was  488 
cells  axially  and  30  cells  radially.  The  computational  and 
experimental  diaphragm  and  test  station  locations  were  the 
same  as  shown  in  Figure  2. 

Five  computations  were  performed  for  the  LBS  configuration 
with  initial  conditions  that  duplicated  experimental  runs. 

The  initial  conditions  are  summarized  in  the  following  table: 


TABLE  1.  1/57  Scale  Model  LBS  Initial  Conditions 


Shot  #  P4/P1 


16.0 


Heated  Driver  T4/T1  Overpressure 

kPa  Cpsi) 

no  1.0  34  (  5) 


12 

136.0 

no 

1.0 

175 

(25) 

85-11 

68.6 

yes 

1.58 

175 

(25) 

85-21 

224.4 

no 

1.0 

241 

(35) 

85-23 

132.1 

yes 

1.84 

241 

(35) 

5.  RESULTS  AND  DISCUSSION 

Computational  results  are  presented  in  this  section  for 
the  five  initial  conditions  listed  in  Table  1.  The  five 
conditions  represent  possible  lower  and  upper  limits  of  blast 
simulation  under  consideration  for  the  LB/TS.  The  175  kPa  and 
241  kPa  pressure  levels  are  more  difficult  to  computationally 
simulate  than  the  34  kPa  level  because  of  the  increased 
complexity  of  the  resulting  flow  gradients.  Comparisons  with 
experimental  data  are  made  where  possible. 


Data  is  presented  in  the  form  of  overpressure  and  dynamic 
pressure  versus  time  plots  at  station  7.  Station  7  is  located 
seven  diameters  downstream  from  the  end  of  the  diverging 
nozzle  and  is  the  primary  test  location.  Also,  density  and 
pressure  contour  plots  are  presented.  The  density  and  pressure 
contour  plots  reveal  the  nature  of  the  flow  physics  in  the 
shock  tube.  In  this  study,  efforts  are  concentrated  on 
understanding  two  flow  phenomena;  one  is  the  nozzle  flow  that 
results  in  the  large  blast  thermal  simulator,  the  second  is 
the  flow  differences  that  result  when  the  driver  gas  is  heated 
versus  when  the  driver  gas  is  not  heated. 

The  pressure  versus  time  plots  for  the  34  kPa  (5  psi) 
run  are  presented  in  Figure  5.  An  experimental  dynamic 
pressure  record  was  not  available  for  this  case,  therefore, 
the  stagnation  pressure  record  was  substituted.  The 
experimental  probe  was  positioned  at  one  half  the  radius  of 
the  tube.  The  computational  data  sampled  at  one  half  the 
radius  reveals  excellent  agreement  while  the  computational 
data  sampled  at  the  centerline  is  in  poor  agreement.  This 
implies  that  the  stagnation  pressure  in  the  radial  direction 
is  nonuniform.  The  stagnation  pressure  versus  time  plot  shows 
the  importance  of  computationally  sampling  at  the  same 
spanwise  position  in  the  tube  as  the  experiment . 
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The  statio  overpressure  plot  shows  a  discrepancy 
between  the  experimental  and  computational  data  on  the 
end  of  the  positive  phase  duration  (time  when  static 
overpressure  goes  negative).  The  experimental  and 
computational  data  indicates  the  end  of  the  positive 
phase  duration  at  26  and  18  ms  respectively.  The  volume 
of  the  computational  driver  after  grid  generation  was 
slightly  smaller  than  the  volume  of  the  experimental 
driver,  which  could  cause  this  discrepancy.  Otherwise, 
the  computational  and  experimental  data  on  the  pressure 
versus  time  plots  are  in  excellent  agreement.  Note  the 
computational  data  is  taken  at  1/4  diameter  and  at  the 
centerline.  The  experimental  data  was  taken  by  probes 
flush  with  the  tube  wall.  This  excellent  comparison 
indicates  that  the  static  overpressure  is  radially 
uniform  in  the  tube.  In  this  paper,  the  computational 
and  experimental  comparison  of  static  overpressure 
produced  very  good  agreement  at  all  shock  overpressures. 

The  contour  data  for  the  34  kPa  (5  psi)  run,  Figures 
6  through  8,  reveal  the  formation  of  a  complex  shock 
system  in  the  diverging  nozzle.  The  contour  data  shows 
the  primary  shock,  contact  surface,  and  backward-facing 
shock  developing  in  the  diverging  nozzle.  The 
backward-facing  shock  is  not  a  normal  shock,  but  consists 
of  two  oblique  shocks  which  terminate  in  a  normal  shock 
at  the  center  of  the  tube.  The  intersection  of  the 
oblique  and  normal  shock  produces  a  reflected  or  an 
oblique  transmitted  shock  which  can  also  be  seen  in  the 
contour  plots.  Figure  7.  The  complex  structure  of  the 
backward-facing  shock  has  been  verified  experimentally, 
as  shown  in  a  shadowgraph  photograph  by  Amann(8),  Figure 
9. 

The  primary  shock  moves  quickly  through  the  diverging 
nozzle  and  proceeds  downstream  as  a  normal  shock.  The  contact 
surface  which  separates  the  driver  gas  from  driven  gas  can  be 
distinguished  as  the  clustered  gradient  present  in  the  density 
plot  and  not  present  in  the  pressure  plot.  The  contact 
surface  moves  downstream,  maintaining  a  curved  surface. 

Behind  it,  the  driver  gas  resembles  core  flow.  In  Figure  8, 
the  contact  surface  takes  on  a  Jetting  appearance  at  the 
centerline  boundary  that  may  be  an  artifice  of  the  code. 

Behind  the  contact  surface,  a  complex  system  of  oblique 


shocks,  rotational  motion  and  slip  surfaces  develops  in  the 
core  flow. 

The  oblique  shocks  are  a  result  of  the  backward 
facing  shock,  which  remains  in  the  diverging  nozzle, 
reflecting  at  the  solid  walls  and  symmetry  boundary  as 
oblique  shocks.  These  reflections  set  up  a  shock  diamond 
pattern  that  stretches  downstream.  This  shock  diamond 
pattern  was  discussed  by  Gottlieb(9),  "If  the  conical 
expansion  is  too  rapid,  however,  or  boundary-layer 
effects  are  significant  this  upstream-facing  shock  wave 
would  be  a  series  of  cells  consisting  of  oblique  shock 
waves.  The  spatial  extent  or  length  of  such  pseudo-shock 
patterns  is  very  large  compared  to  the  thickness  of  a 
typical  normal  shock  wave,  and  they  can  cover  a  duct 
length  of  many  tube  diameters".  The  computations 
presented  here  were  lnviscid,  which  implies  the  rapid 
expansion  alone  can  produce  the  shock  diamond  pattern 
shown  in  the  contour  plots. 

Even  though  this  is  an  lnviscid  code,  rotational  motion  can 
occur,  as  visible  in  Figure  8.  Gradients  in  total  enthalpy 
are  caused  by  the  unsteady  nature  of  the  shock.  These 
gradients  in  entropy  occur  when  some  streamlines  experience  a 
higher  entropy  increase  by  going  through  the  normal  part  of 
the  recompression  shock  while  other  streamlines  experience  a 
lower  entropy  increase  by  going  through  the  oblique  part  of 
the  recompression  shock.  From  Crocco's  theorem  we  know  that 
whenever  gradients  in  total  entalpy  or  gradients  in  entropy 
exist  in  the  flow  field,  rotational  motion  occurs. 

Figures  10  through  13  present  pressure  versus  time 
results  for  the  17b  kPa  (25  psi)  and  241  kPa  (35  psi)  unheated 
and  heated  runs.  As  stated  before,  the  computational  and 
experimental  comparison  of  static  pressure  produced  good 
agreement  at  all  shock  overpressures.  However,  it  is  very 
important  for  blast  simulation  to  model  not  only  static 
pressure  but  dynamic  pressures  as  well. 

Experimental  dynamic  pressure  is  not  measured  directly 
in  the  shock  tube  but  is  calculated  from  static  and  stagnation 
pressure  records.  The  static  pressure  data  was  measured  at 
the  wall  and  the  stagnation  pressure  data  was  measured  at  1/4 
diameter.  To  be  rigorous,  the  static  and  stagnation  probes 


should  15©  at  the  same  location  in  order  to  calculate  true 
local  Mach  number  and  thus  true  dynamic  pressure.  However,  if 
the  static  pressure  is  radially  uniform,  then  there  is  no  need 
to  have  the  probes  at  the  same  exact  radial  location.  The 
code  indicates  static  pressure  is  radially  uniform,  although 
this  has  not  been  confirmed  experimentally  by  placing  probes 
at  different  radial  locations. 

For  the  176  kPa  heated  driver  oase,  the  experimental  and 
computational  comparison  of  dynamio  pressure  are  in  excellent 
agreement  until  the  contact  surface  goes  by  at  .008  seoonds. 
After  the  oontact  surface,  the  computation  at  1/4  diameter 
and  the  experimental  record,  based  on  a  stagnation  probe  at 
1/4  diameter  are  very  similar  except  for  one  anomaly.  The 
motivation  for  heating  the  driver  gas  can  be  seen  in  the 
experimental  records  for  the  175  kPa  unheated  driver  case. 

For  the  175  kPa  unheated  driver  case ,  a  much  noisier 
experimental  dynamic  pressure  record  is  indicated  after  the 
arrival  of  the  contact  surface.  In  the  previous  case,  heating 
smoothed  this  region  and  produced  a  simulation  closer  to  that 
of  a  decaying  free  field  wave.  Computationally,  the  comparison 
is  good  until  the  arrival  of  the  contact  surface.  After  the 
contact  surface,  the  experiment  shows  a  higher  Mach  number 
flow  at  the  centerline. 

The  purpose  of  heating  is  to  reduce  the  jump  in  dynamic 
pressure  (Figures  10  and  12)  which  occurs  across  the  contact 
surface  for  unheated  cases.  The  jump  in  dynamic  pressure 
Increases  with  the  shock  overpressure . (2)  When  the  diaphragm 
is  opened,  the  driver  gas  is  cooled  by  the  passage  of  the 
rarefaction  wave  into  the  driver  section.  The  driven  gas  is 
heated  by  the  passage  of  the  primary  shock.  At  the  contact 
surface,  where  the  driver  gas  meets  the  driven  gas,  the 
difference  in  temperature  of  the  driver  and  driven  gas  results 
in  a  difference  in  density  which  in  turn  is  demonstrated  as  a 
jump  in  the  dynamic  pressure  plots.  By  heating  the  driver 
gas  to  the  proper  level,  the  temperature  on  each  side  of  the 
contact  surface  can  be  matched. 

Furthermore,  heating  of  the  driver  gas  reduces  the  driver 
pressure  required  to  obtain  a  given  shock  overpressure  as 
shown  by  the  initial  conditions  in  Table  1.  Heating  also 
significantly  reduces  the  spike  activity  in  the  stagnation  and 
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dynamic  pressure  versus  time  records  (Figures  11  and  13)  and 
produces  overpressure  versus  time  records  that  more  closely 
resemble  a  smoothly  decaying  blast  overpressure  wave  shape. 

Figures  14  through  19  present  contour  plots  for  the  175 
kPa  (25  psi)  and  241  kPa  (35  psi)  unheated  and  heated  runs. 
Another  effect  of  heating,  as  shown  in  the  last  contour  plots 
(Figures  16  and  19)  available  for  the  time  simulated,  is  the 
backward-facing  shock  is  eventually  swallowed  in  the  nozzle 
for  the  heated  case,  but  remains  outside  the  nozzle  for  the 
unheated  case. 


VI .  Conclusions 

One-dimensional  calculations (3 , 4)  have  had  limited 
success  for  accurately  predicting  the  flow  through  the 
rapid  expansion  of  the  LB/TS  design  because  the  flow  in 
this  region  is  multi-dimensional.  Two-dimensional 
calculations  with  the  upwind,  TVD,  finite-volume, 
implicit  scheme  in  the  BRL  BLAST2D  code  were  presented 
here  which  captured  and  revealed  the  nature  of  the  flow 
physics  in  the  1/57  scale  LBS.  A  complex  system  of 
shocks,  vortices,  and  slip  surfaces  were  revealed  in 
contour  plots. 

Nozzle  flow  produces  a  complex  recompression  shock 
system  which  influences  the  flow  behind  the  contact 
surface.  The  static  pressure  is  uniform  radially, 
however,  the  stagnation  pressure  and  thus  both  the  Mach 
number  and  dynamic  pressure  vary  greatly.  Because  of 
this  radial  variation  in  flow  it  is  very  important  to 
compare  dynamic  pressure  computational  and  experimental 
records  at  the  same  radial  location.  Heating  reduces  the 
required  driver  pressure  for  a  given  shock  overpressure 
and  smooths  the  flow  behind  the  contact  surface, 
producing  a  dynamic  pressure  record  closer  to  that  of  a 
free  field  wave. 

The  BLAST2D  code  provides  excellent  modeling  capability 
at  low  shock  overpressures  and  good  modeling  capabilities  at 
higher  shock  overpressures.  The  code  simulates  heated  driver 
cases  better  than  unheated  driver  cases  for  high  shock 
overpressures.  Future  efforts  will  be  concentrated  on 
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Figure  10:  Pressure  versus  Time  Plots,  175  kPa  (25  psi) 
Peak  Overpressure,  Unheated  Driver,  Shot  #12 
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Figure  15:  Contour  Plots  at  1.5  ms  and  2.0  ms, 
241  kPa  (35  psi)  Peak  Overpressure, 
Shot  #85-21,  Unseated  Driver 


35 


J-V  V'/W  ■' 


v.v; /./.W\V^V'J' 


density  -4  ms 


c 


IM 


pressure  .4  ms 


density  . 7  ms 


c 


w 


pressure  .7  ms 


Figure  18:  Contour  Plots  at  1.1  ms  and  1.5  ms, 
241  kPa  (35  psi)  Peak  Overpressure, 
Shot  #85-23,  Heated  Driver 


Figure  19: 


Contour  Plots  at  4.1  ms  and  17.3  ms, 
241  kPa  (35  psi)  Peak  Overpressure, 
Shot  #85-23,  Heated  Driver 
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